Bouw LHM topologie

Issue: Geen directe koppeling gevonden tussen LSW's en DM via LKM. Daarom hebben we gekozen voor de volgende methode.¶

Methode om LSW's te linken aan DM:¶

1. We vinden LSW eindpunten door routing van LSM links. Aanname:

   • We gaan ervan uit dat LSWs die afwateren op andere LSWs níet tevens afwateren op het DM.

   Resultaat: 2309 LSW eindpunten moeten gekoppeld worden aan DM knopen

2. We checken welke LSW eindpunten gekoppeld zijn aan slechts 1 DM knoop en koppelen deze direct aan de bijbehorende DM-knoop via dwkeys. Methode:

   • Check voor iedere LSW eindpunt het bijbehorende district
   • Check welke districten aan slechts 1 DM-knoop worden gekoppeld door uitlezen van dwkeys.txt

   Resultaat: 1784 LSWs gekoppeld aan DM-knopen

3. De niet gekoppelde LSW eindpunten worden door LSM-lateralen aan DM gekoppeld: Data:

   • We hebben een routing LKM file (LKM25_links) die de afvoerrichting naar LKM netwerk weergeeft (blauwe pijlen in de figuur)
   • We hebben een LKM waterlichamen file (lkm_waterlichamen_shp) (rode polygonen in de figuur)
   • We hebben een LSM3 lateraal file (LSM3_DMKnoopDistrict_childs.csv) die aangeeft welke LSM3 lateral aan DM knoop is gekoppeld (groene punten in de figuur)

   Methode:

   1. We volgen de routing in de LKM file van de LSW eindpunten tot deze uitkomt in een LKM waterlichaam.
   2. We zoeken de dichtsbijzijnde LSM3 lateraal van het LKM waterlichaam
   3. We lezen de bijbehorende DM knoop uit de LSM3_DMKnoopDistrict_childs.csv

   Resultaat: 457 LSWs koppelingen gemaakt aan DM-knopen

   

4. De LSW eindpunten die overblijven (193) koppelen we aan de dichtsbijzijnde DM-knoop

In [1]:

from config import DATA_DIR, LKM25_DIR, MOZART_DIR, LHM_DIR, LSW_DIR, LSM_KOPPELING_DIR, load_src
import geopandas as gpd
load_src()

from lhm.read import read_lsm_lhm, read_lsw_routing, read_dw_keys
from lhm.lsm import snap_to_waterbodies
from lhm.lsw import lsw_end_nodes, lsw_network

from config import DATA_DIR, LKM25_DIR, MOZART_DIR, LHM_DIR, LSW_DIR, LSM_KOPPELING_DIR, load_src import geopandas as gpd load_src() from lhm.read import read_lsm_lhm, read_lsw_routing, read_dw_keys from lhm.lsm import snap_to_waterbodies from lhm.lsw import lsw_end_nodes, lsw_network

1. We vinden LSW eindpunten door routing van LSM links¶

Vanuit de geleverde lsws.shp en lswrouting.dik maken we een netwerk van LSWs.

In [2]:

lsw_routing_dik = MOZART_DIR / r"mozartin/lswrouting.dik"
lsw_routing_df = read_lsw_routing(lsw_routing_dik)
lsw_gdf = gpd.read_file(LSW_DIR / "lsws.shp")

lsw_gdf = lsw_gdf.dissolve(by="LSWFINAL").reset_index()
lsw_links_gdf, lsw_nodes_gdf = lsw_network(lsw_gdf, lsw_routing_df)

lsw_routing_dik = MOZART_DIR / r"mozartin/lswrouting.dik" lsw_routing_df = read_lsw_routing(lsw_routing_dik) lsw_gdf = gpd.read_file(LSW_DIR / "lsws.shp") lsw_gdf = lsw_gdf.dissolve(by="LSWFINAL").reset_index() lsw_links_gdf, lsw_nodes_gdf = lsw_network(lsw_gdf, lsw_routing_df)

Vinden niet-gekoppelde LSWs: Eindpunt LSW's¶

We gaan ervan uit dat LSWs die afwateren op andere LSWs níet tevens afwateren op het DM. We gaan vanaf dit punt LSWs die niet gekoppeld zijn aan andere LSWs koppelen aan het DM. Hieronder tonen we deze LSWs op de kaart

In [3]:

lsw_end_nodes_gdf = lsw_end_nodes(lsw_links_gdf, lsw_gdf)
print(f"we moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()

lsw_end_nodes_gdf = lsw_end_nodes(lsw_links_gdf, lsw_gdf) print(f"we moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen") lsw_end_nodes_gdf.explore()

we moeten 2309 koppelen aan DM-knopen

Out[3]:

Make this Notebook Trusted to load map: File -> Trust Notebook

2. We checken welke LSW eindpunten gekoppeld zijn aan slechts 1 DM knoop en koppelen deze direct aan de bijbehorende DM-knoopKoppelen via DW-keys¶

Vanuit dwkeys.txtkunnen we per district kijken aan welke DM knoop wordt gekoppeld. Wanneer een district slechts koppelt aan 1 DM-knoop, dan betekent dit dat alle LSWs in dit district gekoppeld moeten worden aan deze knoop.

We doen dit en tonen de overgebleven nog niet gekoppelde LSWs.

In [4]:

dw_key_file = LHM_DIR / r"dm/txtfiles_git/dwkeys.txt"

dw_keys_df = read_dw_keys(dw_key_file)
lsw_dm_links = []

for dw, df in dw_keys_df.groupby(by=["oid"]):
    dm_nodes = df.loc[df.kty == "d"].nid.unique()
    if len(dm_nodes) == 1:
        df = lsw_end_nodes_gdf[lsw_end_nodes_gdf["DWRN"] == dw[0]]
        for row in df.itertuples():
            lsw_dm_links += [(row.LSWFINAL, dm_nodes[0])]

print(f"{len(lsw_dm_links)} LSWs gekoppeld aan DM-knopen")
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in lsw_dm_links])]
print(f"We moeten nog {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()

dw_key_file = LHM_DIR / r"dm/txtfiles_git/dwkeys.txt" dw_keys_df = read_dw_keys(dw_key_file) lsw_dm_links = [] for dw, df in dw_keys_df.groupby(by=["oid"]): dm_nodes = df.loc[df.kty == "d"].nid.unique() if len(dm_nodes) == 1: df = lsw_end_nodes_gdf[lsw_end_nodes_gdf["DWRN"] == dw[0]] for row in df.itertuples(): lsw_dm_links += [(row.LSWFINAL, dm_nodes[0])] print(f"{len(lsw_dm_links)} LSWs gekoppeld aan DM-knopen") lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in lsw_dm_links])] print(f"We moeten nog {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen") lsw_end_nodes_gdf.explore()

1784 LSWs gekoppeld aan DM-knopen
We moeten nog 525 koppelen aan DM-knopen

Out[4]:

Make this Notebook Trusted to load map: File -> Trust Notebook

3. De niet gekoppelde LSW eindpunten worden door LSM-lateralen aan DM gekoppeld:¶

Koppelen LSM lateralen aan het DM¶

Vanuit de bestanden LSM3_locations.csv, LSM3_DMKnoopDistrict_childs.csv kunnen we uitzoeken hoe aan welke DM-knoop een LSM-lateraal is gekoppeld. Een LSM-lateraal ligt meestal vrij netjes in het watersysteem geschematiseerd.

In [5]:

lsm3_locations_csv = LSM_KOPPELING_DIR / "LSM3_locations.csv"
knoop_district_csv = LSM_KOPPELING_DIR / "LSM3_DMKnoopDistrict_childs.csv"
lsm_lhm_gdf = read_lsm_lhm(lsm3_locations_csv, knoop_district_csv)

lsm_lhm_gdf.explore()

lsm3_locations_csv = LSM_KOPPELING_DIR / "LSM3_locations.csv" knoop_district_csv = LSM_KOPPELING_DIR / "LSM3_DMKnoopDistrict_childs.csv" lsm_lhm_gdf = read_lsm_lhm(lsm3_locations_csv, knoop_district_csv) lsm_lhm_gdf.explore()

Out[5]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Koppelen LSW via LSM-lateralen aan DM¶

Voor de nog overgebleven LSWs zoeken we de DM-knoop via LSM-lateralen. We nemen deze stappen:

1. We kijken of er LSM-lateralen van het juist district liggen in het gebied van de LSW
2. We kijken aan welke unieke DM-knopen deze lateralen zijn verbonden. We verbinden de LSW met al deze districten

In [6]:

lkm_waterlichamen_shp = LKM25_DIR / "KRW-waterlichamen_SGBP3.shp"
lkm_waterlichamen_gdf = gpd.read_file(lkm_waterlichamen_shp)
lsm_lhm_snapped_gdf = snap_to_waterbodies(lsm_lhm_gdf, lkm_waterlichamen_gdf, offset=250)

lkm_links_shp = LKM25_DIR / "LKM25_Links.shp"
lkm_links_gdf = gpd.read_file(lkm_links_shp)

lkm_waterlichamen_shp = LKM25_DIR / "KRW-waterlichamen_SGBP3.shp" lkm_waterlichamen_gdf = gpd.read_file(lkm_waterlichamen_shp) lsm_lhm_snapped_gdf = snap_to_waterbodies(lsm_lhm_gdf, lkm_waterlichamen_gdf, offset=250) lkm_links_shp = LKM25_DIR / "LKM25_Links.shp" lkm_links_gdf = gpd.read_file(lkm_links_shp)

ToDo: Deze code netjes verwerken in een functie

In [8]:

from shapely.geometry import Point
from lhm.dm import get_dm_nodes

if lsw_gdf.index.name !="LSWFINAL":
    lsw_gdf.set_index("LSWFINAL", inplace=True)
init = len(lsw_dm_links)

def report_progress(iteration, total, interval=10):
    percent = (iteration / total) * 100
    progress = int(percent / interval)
    line = '[' + '#' * progress + ' ' * (int(100 / interval) - progress) + ']'
    print(f'{line} {percent:.1f}% completed', end='\r')

def find_links(row, lsw_id): # algemene functie om links te vinden mbv een LKM link
    lsw_dm_links = [] 
    polygon = lsw_gdf.at[lsw_id, "geometry"]
    df = lsm_lhm_gdf[lsm_lhm_gdf.within(polygon)]# we kijken of er lateralen in het gebied van de LSW liggen
   
    if not df.empty:
        lsw_dm_links +=[(lsw_id, i) for i in df.DM.unique()]
   
    if row.nodeto.startswith("LSM"): # de link moet starten met LSM
        point = Point(row.geometry.bounds[2:]) # de Point van de nodeto
        waterlichaam_gdf = lkm_waterlichamen_gdf[lkm_waterlichamen_gdf.intersects(point)] # het waterlichaam waar punt op snapt
       
        if not waterlichaam_gdf.empty: # er moet wél een waterlichaam gevonden worden
            waterlichaam = waterlichaam_gdf.iloc[0] # dan pakken we het eerste waterlichaam (als het goed is is de lengte altijd 1)
            df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)] # we zoeken naar lsm_lhm laterale knopen binnen het waterlichaam
            df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)]
            
            if not df.empty: # we hebben gevonden!
                lsw_dm_links = [(lsw_id, i) for i in df.DM.unique()] # hier maken we links van de lsw_id naar de unieke DM-knopen   
                lsw_dm_links = [i for i in lsw_dm_links if i[1] in dm_nodes_stringified]
    
    return lsw_dm_links
all_links = []

for idx, row in enumerate(lsw_end_nodes_gdf.itertuples()):
    report_progress(idx, len(lsw_end_nodes_gdf))
    lsw_id = row.LSWFINAL
    district = row.DWRN
    dm_nodes = get_dm_nodes(dw_keys_df,district)
    dm_nodes_stringified = [str(i) for i in dm_nodes]
    lkm_links_iter = lkm_links_gdf[lkm_links_gdf.nodefrom == str(lsw_id)].itertuples() 
   
    for row in lkm_links_iter:
        links = []
        links = find_links(row, lsw_id) 
      
        if (not links):
            counter = 0
         
            while (not links) and (counter < 10) :
                df = lkm_links_gdf[lkm_links_gdf.nodefrom == row.nodeto]
             
                if not df.empty:
                    row = df.iloc[0]
                    links = find_links(row, lsw_id)
                counter += 1
        
        all_links += [i for i in links if i not in all_links]

all_links
lsw_dm_links += all_links
print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen") 
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])]
print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()

from shapely.geometry import Point from lhm.dm import get_dm_nodes if lsw_gdf.index.name !="LSWFINAL": lsw_gdf.set_index("LSWFINAL", inplace=True) init = len(lsw_dm_links) def report_progress(iteration, total, interval=10): percent = (iteration / total) * 100 progress = int(percent / interval) line = '[' + '#' * progress + ' ' * (int(100 / interval) - progress) + ']' print(f'{line} {percent:.1f}% completed', end='\r') def find_links(row, lsw_id): # algemene functie om links te vinden mbv een LKM link lsw_dm_links = [] polygon = lsw_gdf.at[lsw_id, "geometry"] df = lsm_lhm_gdf[lsm_lhm_gdf.within(polygon)]# we kijken of er lateralen in het gebied van de LSW liggen if not df.empty: lsw_dm_links +=[(lsw_id, i) for i in df.DM.unique()] if row.nodeto.startswith("LSM"): # de link moet starten met LSM point = Point(row.geometry.bounds[2:]) # de Point van de nodeto waterlichaam_gdf = lkm_waterlichamen_gdf[lkm_waterlichamen_gdf.intersects(point)] # het waterlichaam waar punt op snapt if not waterlichaam_gdf.empty: # er moet wél een waterlichaam gevonden worden waterlichaam = waterlichaam_gdf.iloc[0] # dan pakken we het eerste waterlichaam (als het goed is is de lengte altijd 1) df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)] # we zoeken naar lsm_lhm laterale knopen binnen het waterlichaam df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)] if not df.empty: # we hebben gevonden! lsw_dm_links = [(lsw_id, i) for i in df.DM.unique()] # hier maken we links van de lsw_id naar de unieke DM-knopen lsw_dm_links = [i for i in lsw_dm_links if i[1] in dm_nodes_stringified] return lsw_dm_links all_links = [] for idx, row in enumerate(lsw_end_nodes_gdf.itertuples()): report_progress(idx, len(lsw_end_nodes_gdf)) lsw_id = row.LSWFINAL district = row.DWRN dm_nodes = get_dm_nodes(dw_keys_df,district) dm_nodes_stringified = [str(i) for i in dm_nodes] lkm_links_iter = lkm_links_gdf[lkm_links_gdf.nodefrom == str(lsw_id)].itertuples() for row in lkm_links_iter: links = [] links = find_links(row, lsw_id) if (not links): counter = 0 while (not links) and (counter < 10) : df = lkm_links_gdf[lkm_links_gdf.nodefrom == row.nodeto] if not df.empty: row = df.iloc[0] links = find_links(row, lsw_id) counter += 1 all_links += [i for i in links if i not in all_links] all_links lsw_dm_links += all_links print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen") lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])] print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen") lsw_end_nodes_gdf.explore()

457 LSWs koppelingen gemaakt aan DM-knopen
We moeten 193 koppelen aan DM-knopen

Out[8]:

Make this Notebook Trusted to load map: File -> Trust Notebook

4. Koppelen LSW aan dichtsbijzijnde DM-knoop¶

Als uiterste mogelijkheid koppelen we LSW-knopen aan de meest dichtbijzijnde DM-knoop die voorkomt in de DM-keys

Locatie DM knopen inlezen

In [9]:

# Define the shapefile paths
DM_nodes_shp = LHM_DIR / r"dm/data_dvc/DM_nodes.shp"
DM_links_shp = LHM_DIR / r"dm/data_dvc/DM_links.shp"

# Read each shapefile separately
DM_nodes_gdf = gpd.read_file(DM_nodes_shp)
DM_links_gdf = gpd.read_file(DM_links_shp)
DM_nodes_gdf.columns = [i.strip() for i in DM_nodes_gdf.columns]
DM_links_gdf.columns = [i.strip() for i in DM_links_gdf.columns]

# Calculate the closest distance LSW endpoints and DM
all_links = []
for row in lsw_end_nodes_gdf.itertuples():
    lsw_id = row.LSWFINAL
    district = row.DWRN
    dm_nodes = get_dm_nodes(dw_keys_df,district)
    dm_nodes = [str(i) for i in dm_nodes]

    dm_node = DM_nodes_gdf.at[DM_nodes_gdf.loc[DM_nodes_gdf.ID.isin(dm_nodes)].distance(row.geometry).sort_values().index[0], "ID"]

    all_links += [(lsw_id, int(dm_node))]

lsw_dm_links += all_links
print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen") 
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])]
print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")

# Define the shapefile paths DM_nodes_shp = LHM_DIR / r"dm/data_dvc/DM_nodes.shp" DM_links_shp = LHM_DIR / r"dm/data_dvc/DM_links.shp" # Read each shapefile separately DM_nodes_gdf = gpd.read_file(DM_nodes_shp) DM_links_gdf = gpd.read_file(DM_links_shp) DM_nodes_gdf.columns = [i.strip() for i in DM_nodes_gdf.columns] DM_links_gdf.columns = [i.strip() for i in DM_links_gdf.columns] # Calculate the closest distance LSW endpoints and DM all_links = [] for row in lsw_end_nodes_gdf.itertuples(): lsw_id = row.LSWFINAL district = row.DWRN dm_nodes = get_dm_nodes(dw_keys_df,district) dm_nodes = [str(i) for i in dm_nodes] dm_node = DM_nodes_gdf.at[DM_nodes_gdf.loc[DM_nodes_gdf.ID.isin(dm_nodes)].distance(row.geometry).sort_values().index[0], "ID"] all_links += [(lsw_id, int(dm_node))] lsw_dm_links += all_links print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen") lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])] print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")

193 LSWs koppelingen gemaakt aan DM-knopen
We moeten 0 koppelen aan DM-knopen

In [10]:

from shapely.geometry import LineString
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
lsw_dm_links_gdf = gpd.GeoDataFrame(lsw_dm_links, columns = ["node_from", "node_to"], geometry = gpd.GeoSeries(), crs=28992)

def make_line_string(row):
    report_progress(row._name, len(lsw_dm_links_gdf))
    point_from = lsw_nodes_gdf.loc[row.node_from]
    point_to = DM_nodes_gdf.loc[DM_nodes_gdf.ID == str(row.node_to)].geometry
    return LineString([[point_from.x, point_from.y], [point_to.x, point_to.y]])

lsw_dm_links_gdf.loc[:, "geometry"] = lsw_dm_links_gdf.apply((lambda x: make_line_string(x)), axis=1)

lsw_dm_links_gdf["node_from"] = lsw_dm_links_gdf["node_from"].astype(str)
lsw_dm_links_gdf["node_to"] = lsw_dm_links_gdf["node_to"].astype(str)

from shapely.geometry import LineString import warnings warnings.simplefilter(action='ignore', category=FutureWarning) lsw_dm_links_gdf = gpd.GeoDataFrame(lsw_dm_links, columns = ["node_from", "node_to"], geometry = gpd.GeoSeries(), crs=28992) def make_line_string(row): report_progress(row._name, len(lsw_dm_links_gdf)) point_from = lsw_nodes_gdf.loc[row.node_from] point_to = DM_nodes_gdf.loc[DM_nodes_gdf.ID == str(row.node_to)].geometry return LineString([[point_from.x, point_from.y], [point_to.x, point_to.y]]) lsw_dm_links_gdf.loc[:, "geometry"] = lsw_dm_links_gdf.apply((lambda x: make_line_string(x)), axis=1) lsw_dm_links_gdf["node_from"] = lsw_dm_links_gdf["node_from"].astype(str) lsw_dm_links_gdf["node_to"] = lsw_dm_links_gdf["node_to"].astype(str)

[######### ] 100.0% completed

Resultaat¶

In de kaart hieronder zie je het resultaat

In [11]:

import folium
#m = lsw_gdf.explore(name="lsws", color="lightgrey")
m = lsw_links_gdf.explore(name="lsw links", color="blue")
m = lsw_dm_links_gdf.explore(m=m, name="lsw-dm links", color="green")
m = DM_links_gdf.explore(m=m, name="dm links", color="red")

folium.LayerControl().add_to(m)

m

import folium #m = lsw_gdf.explore(name="lsws", color="lightgrey") m = lsw_links_gdf.explore(name="lsw links", color="blue") m = lsw_dm_links_gdf.explore(m=m, name="lsw-dm links", color="green") m = DM_links_gdf.explore(m=m, name="dm links", color="red") folium.LayerControl().add_to(m) m

Out[11]:

Make this Notebook Trusted to load map: File -> Trust Notebook